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ABSTRACT 

This paper presents RK method based on Heronian mean for solving Delay differential equations with constant 
lags (delays). The delay term is approximated by using linear and Lagrange interpolation. Two numerical examples have 
been considered to illustrate the effectiveness of this method. The numerical results are also compared with RK method 
based on Arithmetic mean. 
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1. INTRODUCTION 

Delay differential equations (DDEs) arise in chemical kinetics [1], control systems [2], population dynamics [3] 
and in many areas of Science and engineering. Recently there has been great interest in finding the numerical solutions of 
DDEs. Many numerical methods have been suggested for solving DDEs of first order as in [4-6], Paul and Baker [7] 
discussed about the determination of stability region of RK method for DDEs. Hu et al. [8] considered the stability of RK 
methods for DDEs with multiple delays. 

Evans and Yaakub [9] have proposed the fourth order RK method based on Heronian mean (RKHeM) for solving 
IVPs. Ponalagusamy and Senthilkumar [10] explained a new RK embedded method with error control. 

In this paper the fourth order RKHeM method has been discussed for obtaining the numerical solution of first 
order DDEs. Here we consider the first order DDEs in the following form: 

y V) = t), y(t -t)), t>t 0 

y(t) = 0(t), t<t 0 

If it has one delay only, or 

yV) = t > t 0 

y(t) = Q>(t), t<t 0 

If it has more than one delay term where O(T) is the initial function? Here the delay terms r, T t , r 2 ,.„, T n arc 
positive constants. Two numerical examples have been considered to demonstrate the adaptability of RKHeM to DDEs. 
This paper has been organized as follows: 
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In Section 2, the fourth order RKAM and RKHeM formula for solving ODEs have been briefly mentioned. In 
Section 3, the adaptability of RKHeM formula for solving DDEs has been discussed. In Section 4, the convergence and 
stability for RKHeM have been discussed. In Section 5, two numerical examples have been provided to illustrate the 
effectiveness of RKHeM. 


2. RKHeM METHOD FOR ORDINARY DIFFERENTIAL EQUATIONS 

Consider the first order equation of the form 

y' = f(x,v) With v(x 0 ) = y 0 . 

Hie fourth order RKAM formula is given by 


y „ + 1 = y„ + - 


h 


k 3 + k 2 + k 2 + k 3 ^ k 3 + k 4 


3 i=l ^ $ i= i 


Where 


K =f{x n ,y n ) 


k 2 = f\ x n + y„ + 2 hk x 


k i = f \ X n y„ + \ hk 2 


k 4 = f( x „ + h, y n + hk 3 ) 

Hie Heronian means of two points y l and y 2 is defined as 
HeM = + y + ^/xv J . 

By replacing the AM in eqn. (2. 1) by HeM, the fourth order RKHeM formula is written as 


K =f(x n ,y„) 
k 2=f( x „ +\>yy„+\h k ,) 

1 1 25 

k i = /(*„ + - ■ K y„ - — hk x + — hk 2 ) 

Z 4o 4o 

1 47 289 

k 4 = f(x n +h,y hk , + hk 2 + hk 3 ) 

4 " " 24 600 300 3 

y „ + 1 = y n +^h[k 1 +2(k 2 +k 3 ) + k 4 +J^ + J^+Jk^~ (or) 


y „+ 1 - y„ +~h k t +2(k 2 +k 3 ) + k 4 + 3 J\k l k 2 \ + ^\ k 2 k 3\ + yj\ k 3 k 4 
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To avoid the negative values under the square root signs. 

3. RKHeM METHOD FOR DELAY DIFFERENTIAL EQUATIONS 

Consider the first order DDE of the form 

y\t) = t>t 0 


y(t) = 0(t), t<t 0 


Where ®(t) the initial is function and T is the delay term. When we adapt the fourth order RKHeM formula to 
DDE, we have 


y n + 1 ~y n + ~^ fl K + 2 (^2 + £3 ) + k A + yj\ k x k 2 1 + yj\ k 2 k 3 1 + J k-,k 4 


where 


w 


+ £,h,y > 


+ h± 

j = i 


a^r 


y(x n +c t h-T) 


(i = 


The above RKHeM formula can also be extended to solve DDEs with multiple delays. 

4. CONVERGENCE AND STABILITY OF RKHeM 

4.1 Order of Convergence of RKHeM 


When we use RKHeM formula to solve DDEs, the delay term y(x n +Cjh — Z") need to be interpolated to 

approximate the value. Many techniques are available for the approximation. In this paper linear and Lagrange 
interpolation are used to approximate the delay term. The interpolation has to be adapted to the order of the method. 

For any given RK method, its adaptation to DDEs by means of interpolation procedure has an order of 
convergence equal to min j p,q} where p denotes the order of consistency of the RK method and q is the number of 
support points of the interpolation procedure. (See [ 11 ]) 

For linear interpolation we are using two support points so that the order of convergence of RK method is two 
(See [ 12 ]) which is less than the order of the RK method. For Lagrange interpolation we are using five support points so 
that the order of convergence of RK method is four for DDEs also. 

4.2 Linear Stability for RKHeM: 

There are many concepts of stability of numerical methods when applied to DDE, depending on the test equation 
as well as the delay term involved. Here our attention is to a linear test equation with a constant delay T , 


y\x) = Ay(x) + juy(x-z), t>t 0 


y(t)=®(t), t<t 0 

Where X, /J £ R , T > 0 and® is continuous. 
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The stability function of s-stage RK method for ODE is given by 


r( z) = 1 + zb r (/ - zA ) 1 e, 


Where 


z = hA, e = (l,l,...,l) r , b = (b l ,b 2 ,...b s f , A = [a..y 

In the case of DDEs, we consider T — Nh where N is a positive integer. In this case when U (^nh + C t h — Z j is 
required a previous internal stage-value Y n _ N ■ is used. Here we consider the stability properties of a recurrence of the 


y n+ 1 = 1 + /L hb T (I - lhA) 1 e y n + /jhb T (I -XhA) l u n _ N 

Where is a vector consisting of ‘back- values’ U ( nh + C f h — r) . 

This can be expressed in a convenient form as 

y n+ 1 = r(Ah)y n + juhb T (l- AhA)~ l u n _ N . 

Ifwe write S = S(Ah) = (/ — AhA) 1 and O n = ^Y n 15 ..., Y n v J , the above can be written as 

y n+ 1 = r Wy n + ^b T S(p n _ N 

We can express this as the recurrence: 

® u =x<b^+z<s> 9 _ N 

Where 

(r(M) (A f o juhb T s) 

A = And Z = — 

1, Se 0) lo juhSA 


This recurrence is stable if the zeros C, [ of the stability polynomial 

S h {X^q)=tei\fl-c?- x X-Z_. 

The root condition for stability is the requirement that all the zeros of S h (/l, /l\C^) satisfy k"/ 1 — 1 , and if 


l 5 ,- 1 then is semi-simple. The stability polynomial of RKAM i 


is obtained as, 
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5 N _5AM 


S(Ah, M h;g)=g 3 -g 


2 6 24 






JSN-\ 




(juhf 


f 


1 + 2 h + 


(Ah) 


2 A 


V 


.2 N- 


(M) 




y 


V 24 y 


And the stability polynomial of RKHeM is obtained as, 

( 


SWi,iHg)=? -g 5N ~ x 


1 +Ah+ 


(Ah)- (Ah) (Ah) mi(Ah) 


5 \ 


-+■ 


24 331776 


J 


-g 4N 1 ( juh) 1+- 2/z+— ( AhY 
v \ 3 96 v y J 


O 25, A 


-g N ~'(vh)- -+-{M) -g 2N -\vh) 

V3 96 y 


^ 25 ^ 

v288 y 


5. NUMERICAL EXAMPLES 
Example 1 

Consider the first order DDE with single delay 

y\t) = yit- 1 ) 

y(t) = Ae f , -l<t<0 


With exact solution is, 


y(t) 


e 

< A + - + e ( - 2+t) A + tA~ — 
e e 

2A- — + e^ +t) A + — - — 
e 2 2e 


on [0,1] 
on [1,2] 
on [2,3] 


This problem is solved by RKAM and RKHeM by using linear interpolation and Lagrange interpolation for the 
delay term with 2 = 1. The absolute error results are shown in the following Tables 1 and 2. 

Table 1: Results of Example 1 (Linear Interpolation) 


Time 

Absolute Error 
(RKAM) 

Absolute Error 
(RKHeM) 

0.50 

1.99e-006 

1.82e-006 

1.00 

5.27e-006 

4.83e-006 

1.50 

7.71e-006 

7.16e-006 

2.00 

1.27e-005 

1.19e-005 

2.50 

1.79e-005 

1.66e-005 

3.00 

2.61e-005 

2.42e-005 
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Table 2: Results of Example 1 (Lagrange Interpolation) 


Time 

Absolute Error 
(RKAM) 

Absolute Error 
(RKHeM) 

0.50 

3.89e-013 

4.12e-012 

1.00 

5.27e-006 

1.21e-010 

1.50 

5.27e-006 

1.17e-010 

2.00 

5.27e-006 

1.08e-010 

2.50 

7.90e-006 

1.64e-010 

3.00 

1.05e-005 

2.17e-010 


Example 2 

Consider the system of first order DDE with multiple delays 

y\(t)=y 5 (t-l)+y 3 (t-l); y , 2 (t)=y l (t-l)+y 2 (t-^); y\(t)=y 3 (t-l)+y l (t-^); 

y\(t)=y 5 (t-l)y 4 (t-l); y' 5 (t)=y l (t-l) for f>0 

With initial functions 

y^(t) = e (l+1) ; y 2 (t) = e ,+l) ; y 3 (0 = sin 0 + 1); 

y 4 (t) = e (t+1) ; y 5 (t) = e {t+l) for t< 0. 

With analytical solutions 


y ] (t)=e f -cost+e; y 2 (t)=2e? +e 2 -2; 
y 3 (t) =e t+2> —cos/ +1 — e 2 +sin(l); 

y 5 (t)=J+e - 1 for 0<f<^. 

This problem is solved by RKAM and RKHeM by using linear interpolation and Lagrange interpolation for the 
delay term. The absolute error results are shown in the following Tables 3 -6. 


Table 3: Results of Example 2 (Linear Interpolation) 


Time 

Absolute Error (RKAM) I 

Yi 

yz 

ys 

yr 

ys 

0.10 

8.35e-007 

1.75e-006 

1.40e-006 

1.85e-006 

8.76e-007 

0.20 

1.68e-006 

3.69e-006 

2.88e-006 

4.10e-006 

1.85e-006 

0.30 

2.54e-006 

5.83e-006 

4.44e-006 

6.85e-006 

2.92e-006 

0.40 

3.44e-006 

8.20e-006 

6.10e-006 

1.02e-006 

4.10e-006 

0.50 

4.39e-006 

1.08e-006 

7.89e-006 

1.43e-006 

5.41e-006 
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Table 4: Results of Example 2 (Linear Interpolation) 


Time 

Absolute Error (RKHeM) 

yi 

Y2 

ya 

yr 

ys 

0.10 

5.70e-007 

1.61e-006 

1.1 le-006 

1.54e-006 

8.03e-007 

0.20 

1.17e-006 

3.38e-006 

2.30e-006 

3.42e-006 

1.69e-006 

0.30 

1.80e-006 

5.35e-006 

3.58e-006 

5.7 le-006 

2.67e-006 

0.40 

2.48e-006 

7.51e-006 

4.97e-006 

8.5 le-006 

3.76e-006 

0.50 

3.21e-006 

9.91e-006 

6.48e-006 

1.19e-005 

4.96e-006 


Table 5: Results of Example 2 (Lagrange Interpolation) 


Time 

Absolute Error (RKAM) 

yi 

ya 

ya 

yr 

ys 

0.10 

2.25e-013 

8.77e-012 

3.07e-012 

9.23e-012 

4.38e-012 

0.20 

9.52e-013 

1.85e-01 1 

6.94e-012 

2.05e-011 

9.23e-012 

0.30 

2.27e-012 

2.92e-01 1 

1 . 1 7e-0 1 1 

3.43e-01 1 

1.46e-01 1 

0.40 

4.28e-012 

4.10e-01 1 

1.76e-011 

5. 1 le-01 1 

2.05e-011 

0.50 

7.07e-012 

2.93e-008 

1.43e-007 

7.16e-01 1 

2.70e-01 1 


Table 6: Results of Example 2 (Lagrange Interpolation) 


Time 

Absolute Error (RKHeM) 

yi 

ya 

ya 

yr 

ys 

0.10 

6.62e-012 

3.64e-012 

7.26e-012 

7.70e-012 

1.82e-12 

0.20 

1.27e-01 1 

7.65e-012 

1.43e-01 1 

1.71 e-0 1 1 

3.82e-012 

0.30 

1.85e-01 1 

1.21e-01 1 

2. 1 3e-0 1 1 

2.85e-011 

6.04e-012 

0.40 

2.39e-01 1 

1.70e-01 1 

2.82e-01 1 

4.25e-01 1 

8.49e-012 

0.50 

2.92e-01 1 

9.56e-01 1 

3.23e-010 

5.96e-011 

1 . 12e-0 1 1 


CONCLUSIONS 

In this paper, the fourth order RKAM and RKHeM formula have been adopted to solve the delay differential 
equations with constant lags. The effectiveness of this approach has been illustrated via examples of DDE with single delay 
and multiple delays. To interpolate the delay term, both linear interpolation and Lagrange interpolation have been 
considered here. 

From the numerical results, it is observed that the RKHeM method is well suitable for solving DDEs. It also 
suggests that the best results can be obtained when we use Lagrange interpolation with suitable number of support points 
for getting fourth order convergence. 
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